% *****************************************************************************
%%                                                                            *
%%                                                                            *
% AUTHOR: Vigneshwar Ramakrishnan											  *	
% Date created: 29 August 2022												  *	
% CC-BY-SA-NC                                                                 *
%Calculation of translation density of states                                 *
%******************************************************************************
% Fourier Transformation
%               N  
% X(k) =       sum  x(n)*exp(-j*2*pi*(k-1)*(n-1)/N), 1 <= k <= N.
%              n=1 
% Constants

clear;clc;
warning off all

kb = 1.3806503E-23;                 % boltzmann constant: m^2 kg s^-2 K^-1
T = 300;                           % Temperature: K
NA = 6.023E23;                     % Avogadro number
mwt_water = 0.018;                 % mol. weight of water: kg/mol (18 g/mol)
m_water = mwt_water/NA;            % mass of 1 water molecule: kg     
beta = 1/(kb*T);                    % used in w_s_HO calculation
h = 6.626068E-34;                  % m^2 kg s^-1
% N_water = 376;
N_water = 128;
%  rho = 20.1242E27 % based on 4 Ang volume thickness (for protein-DNA)
% rho = 21.4106E27; % for water around 3.5 Ang of CAATTC-free DNA
%rho = 91.203E27; % bulk density CAATTC-free DNA
% rho = 22.0968E27 % Density of water around 4 Ang of free CAATTC DNA
 rho = 15E27 % Density of water around 4Ang of protein only.
 %rho = 30.577E27 % Bulk density of water around protein only.
% rho = 30.218E27 % Bulk density of water around protein-DNA complex

% Input Section

a = textread('vac-200-300ps.txt');


t = a(:,1);
x = a(:,2);
n2 = length (x); 

% Modify Sampling

n3 = 2;

m = round(n2/n3);

j = 1;
for i = 1:m

 time(i) = t(j);
 x_mod1(i) = x(j);

 j = j + n3;
 
end

% Make the VAC symmetric

% Negative time

x_symm(1) = x_mod1(m);
for di = 2:m
    x_symm(di)=x_mod1(m-(di-1));
end

% positive time

for dj = 1:m-1
    x_symm(m+dj)=x_mod1(dj+1);
end
    
N = length(x_symm);

% Premultiply by the constants
x_mod = x_symm*N_water*m_water*1E6*(2/(kb*T));

% Declare variables
hist_sum = 0;
trap_sum = 0;

% Frequency Calculations

dt = time(2)-time(1);
Ws = 1/(dt);
freq = Ws*(((-N+1)/2):((N-1)/2))/N;
%freq_sav = (freq((N+1)/2:N)/3E10)';
save('freq-200-300ps.txt','freq','-ascii')
%freq_causal = Ws*(0:N-1) /N;


% **** Non-Causal Fourier Transform Routine

count = 1;

for k = ((-N+1)/2):((N-1)/2)
    
% for k = 1:round(N/2)-1
     
    counter = 1;  
    
    for n = ((-N+1)/2):((N-1)/2)
       
        temp(counter) = x_mod(counter)*exp(-1i*2*pi*(k)*(n)/N);
        counter = counter+1;
        
    end

    % trapezoidal integration of temp
    
    trap_sum = 0;
    
    trap_sum = (temp(1)+temp(N))/2;
    for i = 2:N-1
        trap_sum = trap_sum + temp(i);
    end
           
    X(count) = trap_sum*dt;  
    
    % Histogram integration of temp
    
    hist_sum = 0;
    
    for i = 1:N
        hist_sum = hist_sum+temp(i);
    end
    
    X_hist(count) = hist_sum*dt;
  
count=count+1
end

mod = abs(X);
%mod1 = abs(X_symm);

Fp = mod(1:N);
%Fp1 = mod1(1:N);
save('fp-200-300ps.txt','Fp','-ascii');
%plot (freq,Fp,freq_causal,Fp1,'.');
%fig = plot (freq,Fp,freq,Fp1,'.');
%xlabel ('Frequency, (ps^-1)')
%ylabel ('g , (ps)')
%title ('Density of states for 2000 water')
% 
%saveas(h,'dos.jpg');

% g_t0 calculation

g_t0 = Fp((N+1)/2)*1E-12; % Unit Conversion ps to s

D = g_t0*kb*T/(12*m_water*N_water)


% delta calculation

% rho = 40.62425*6.023*1e23; 


delta = (2*g_t0/(9*N_water))*(((pi*kb*T)/m_water)^0.5)*(rho^(1/3))*((6/pi)^(2/3)) %
% % 
% % % f calculation

fluid = newton(0.3,delta);

% g_g calculation

% 
% delta = 0.0278;
% fluid = 0.3269;

for ci = 1:N
    
    g_g(ci) = Fp((N+1)/2)/(1+((pi*(Fp((N+1)/2)*freq(ci)))/(6*fluid*N_water))^2);
    
end

% g_s calculation

g_s = Fp - g_g; 

%%%%% Entropy calculation %%%%%% 

% Weighting function for Gas-like component

% w_gs = S_HS/(3*k); 

 y = (fluid)^(5/2)/(delta^(3/2)); 

zy = (1+y+(y^2)-(y^3))/(1-y)^3; 

s_hs = ((5/2)+log((((2*pi*m_water*kb*T)/(h^2))^(3/2))*(1/(fluid*rho))*zy)+((y*(3*y-4))/(1-y)^2))*kb % see Lin et al 2010 eq 16

w_gs = s_hs/(3*kb);

% Weighting for solid-like component

bhw = beta * h * freq*1E+12;

for i = ((N+1)/2):N % 
    w_s_HO(i) = (bhw(i)/(exp (bhw(i)) - 1)) - log (1-exp (-bhw(i))); 
end 


% Integrating the weighted density of states

% solid component



for ai = ((N+1)/2):N
    
    solid(ai) = g_s(ai)*1E-12*w_s_HO(ai);

end

solid((N+1)/2) = 0;

freq_ins = freq*1E+12;

P_solid = trapz(freq_ins(((N+1)/2):N),solid(((N+1)/2):N));

for bi = 1:N
    
    gas(bi) = g_g(bi)*1E-12*w_gs;
end

P_gas = trapz(freq_ins(((N+1)/2):N),gas(((N+1)/2):N));

Entropy = P_solid+P_gas; % For N_water water molecules

Finalvalue = Entropy*kb*NA/N_water % J mol^-1 K^-1

% Plot total, solid and gas components

plot(freq/3E-2,Fp*3E-2,'r',freq/3E-2,g_g*3E-2,'g',freq/3E-2,g_s*3E-2,'b')
freq_sav = 
save('freq-200-300ps.txt','Fp','-ascii');
save('fp-200-300ps.txt','Fp','-ascii');

% Energy Calculations

% Weight functions

for ki = 1:N
    w_E_HO(ki) = (bhw(i)/2)+(bhw(i)/(exp(bhw(i))-1));
end

for li = ((N+1)/2):N
    
    E_solid(li) = g_s(ai)*1E-12*w_E_HO(li);

end

E_gas = g_g*1E-12*0.5;

E_gas = trapz(freq_ins(((N+1)/2):N),E_gas(((N+1)/2):N));
E_solid = trapz(freq_ins(((N+1)/2):N),E_solid(((N+1)/2):N));

Energy = E_gas+E_solid;

Final_energy = Energy*NA/N_water*(kb*T);

% w_E_gas = 0.5;

% E = V0 + beta^-1
 


